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Due to the weak spin-orbit interaction and the pecuhar relativistic dispersion in graphene, there 
are exciting proposals to build spin qubits in graphene nanoribbons with armchair boundaries. 
However, the mutual interactions between electrons are neglected in most studies so far and thus 
motivate us to investigate the role of electronic correlations in armchair graphene nanoribbon by 
both analytical and numerical methods. Here we show that the inclusion of mutual repulsions leads 
to drastic changes and the ground state turns ferromagnetic in a range of carrier concentrations. Our 
findings highlight the crucial importance of the electron-electron interaction and its subtle interplay 
with boundary topology in graphene nanoribbons. Furthermore, since the ferromagnetic properties 
sensitively depends on the carrier concentration, it can be manipulated at ease by electric gates. 
The resultant ferromagnetic state with metallic conductivity is not only surprising from an academic 
viewpoint, but also has potential applications in spintronics at nanoscale. 



PACS numbers: 73.22.-f, 72.80.Rj, 75.70.Ak 
I. INTRODUCTION 

Graphene,^ a single hexagonal layer of carbon atoms in 
two dimensions (2D), is the building block for graphitic 
materials ranging from OD fullerenes to ID nanotubes, 
and also the commonly seen 3D graphite. Since it 
was generally believed that the two-dimensional lattice 
should not exist at finite temperature, graphene had of- 
ten been used as a toy model and viewed as an aca- 
demic material until its recent discovery in laboratory.*^ 
The honeycomb structure gives rise to linear dispersion, 
making electrons and holes in graphene massless as in 
relativistic theories.!^ Therefore, most studies focus on 
the electronic and trans port pr operties arisen from this 
peculiar band structurej ^ l ^ l ^ l ^ l ^ l such as the half-integer 
quantum Hall effeciP^ due to the tt Berry phase, the 
quantization of minimal conductivity where carriers are 
almost absentP and so on. 

One of the potential applications of graphene is to 
realize fast electronics at nanoscale, making graphene 
nanoribbons (GNRs) a natural building block for these 
devices. Since the electronic structure sensitively de- 
pends on the transverse width an d also the edg e topology, 
there are intensive investigation^^SEHIlIIlEl] on narrow 
GNRs with width less than 10 nm. Although GNRs have 



been successfully fabricated by lithographji^^down to the 
widths of 20 nm, the roughness of the edges remains large 
(about 5 nm or larger). As a result, theoretical predic- 
tions may not be applicable and limit the fundamental 
and practical applications. A recent breakthrough of 
fabricating GNRs came from chemical methods.!^ It is 
rather remarkable that the width of the GNRs can be 
fabricated in a controlled way down to 10 nm. In ad- 
dition, the edges of these GNRs are ultra smooth with 
possibly well-defined zigzag or armchair shapes, suitable 
for building electric junctions at molecular scale. 

As the transverse width shrinks, the quantum fluctu- 
ations become important and results/predictions from 
mean-field theories shall be checked carefully. Mean- 
while, since the open boundaries of GNR play a crucial 
role at nanoscale, the interplay between the Coulomb 
interaction and the edge morphology will lead to rich 
physics. For instance, it has been revealed that the 
Coulomb interaction gives rise to edge moments in zigzag 
GNR.^^ Furthermore, Son, Cohen, and Louie^*^ have 
shown that, in the presence of external electric field in 
the transverse direction, the system turns half metal- 
lic with magnetic properties controlled by the external 
electric field. Their results not only show that the elec- 
tronic spin degrees of freedom can be manipulated by 



the electric fields, but also bring up the possibility to ex- 
plore spintronics^^ ^■^ ^■^ at the nanometer scale based on 
graphene. 

Inspired by these discoveries, we investigate the effect 
of Coulomb interaction in armchair GNR as schemati- 
cally shown in Fig. [T| Note that the edges are hydro- 
genated so that the dangling a bonds are saturated and 
only the tt bands remain active in low energyP^I By com- 
bining analytical weak-coupling analysis, numerical den- 
sity matrix renormalization-group (DMRG) method, and 
the first-principles calculations, we show that armchair 
GNR exhibits an interesting carrier-mediated ferromag- 
netism upon appropriate doping. Even though only tt 
bands are active in low-energy, in appropriate doping 
regimes, the armchair edges give rise to both itinerant 
Bloch and localized Wannier orbitals. As will be ex- 
plained later, these localized orbitals are direct conse- 
quences of quantum interferences in armchair GNR and 
form flat bands with zero velocity. The carrier-mediated 
ferromagnetism can thus be understood in two steps: 
Electronic correlations in the flat band generate intrinsic 
magnetic moments first, then the itinerant Bloch elec- 
trons mediate ferromagnetic exchange coupling among 
them. As a result, the magnetic properties of armchair 
GNR sensitively depend on the doping and thus can be 
manipulated easily by the external electric fields. 

Though the ferromagnetic state in armchair GNR 
stems from the flat-band states which are partially filled, 
the mechanism is different from the "flat -band fe rromag- 
netism" proposed by Mielke and Tasaki.l^S l^^l^^l The key 
to Mielke- Tasaki ferromagnetism is the finite overlap of 
adjacent Wannier orbitals in the fiat band: the finite 
overlaps generate exchange coupling among these orbitals 
and lead to the flat-band ferromagnetism. However, the 
Wannier orbital in armchair GNR (shown in Fig. [ij has 
zero overlap with its adjacent neighbors. The flat band 
alone only accounts for the existence of the magnetic mo- 
ments and the ferromagnetic order sets in only when itin- 
erant carriers are present. A similar mechanism of ferro- 
magnetism has been arg ued fo r the Hubbard model in a 
kind of frustrated lattice. I^^I^Il 

In the following, we will elaborate in details how the 
carrier-mediated ferromagnetism emerges (upon appro- 
priate doping) in the armchair GNR. In Sec.|ll] we start 
with the Hubbard model and solve the band structure 
by the method of generalized Bloch theorem. In Sec. 
|III| we integrate out the gapped modes and explore the 
ground state properties in weak coupling. We first show 
how the local moments in the fiat band form from the 
Coulomb interaction. We also demonstrate the crucial 
role of itinerant carriers in dispersive bands, which medi- 
ate the indirect exchange coupling among the local mo- 
ments and give rise to the ferromagnetic ground state. 
Indeed, without those carriers, the ferromagnetism disap- 
pears and only Curie-like susceptibility remains in arm- 
chair GNR. In Sec. |IV[ we employ the technique of non- 
Abelian DMRG to investigate the ferromagnetic ground 
state in intermediate coupling. It is remarkable that the 
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FIG. 1: (Color online) (a) Armchair GNR of = 4 and 
Ly = 5. Open edges are present aX y = 1 and Ly. The solid 
circles and squares represent sublattice A and B respectively. 
The shaded circles show the amplitudes of a local Wannier 
orbital oi E = t aX x = 2, with opposite signs indicated by 
light/dark colors. Band structures for the infinitely long GNR 
with (b) Ly = 5 and (c) Ly — 7 clearly show the pair of flat 
bands at E = iit. 

numerical results agree with the weak-coupling analy- 
sis pretty well. In Sec. |Vj the realistic band structure 
and the long-range Coulomb interactions are included via 
the first-principles calculations. It is rather unexpected 
that the fiat band is robust even when the realistic band 
structure is taken into account. Ferromagnetism appears 
around the same doping regime as predicted by either 
weak-coupling or DMRG approaches. Finally, in the last 
section, we discuss the robustness of our predictions and 
also their connections to practical experiments. 

II. HUBBARD MODEL FOR ARMCHAIR 
NANORIBBON 

To understand the carrier-mediated fiat-band ferro- 
magnetism in the armchair GNR, we start with the Hub- 
bard Hamiltonian, 

H = -t J2 [4(r)cJr')+H.c.] + C/^nT(rK(r),(l) 

(r,r'),a r 

where t is the hopping amplitude on the honeycomb net- 
work, t7 > is the on-site repulsion, a =t,i is the spin 
index. The lattice points r = (x,y,A) are labeled by in- 
teger indices {x, y) in longitudinal and transverse direc- 
tions and the additional sublattice index A = A,B. The 
transverse integer index y — 1,2, ...,Ly defines the width 
of the ribbon while x = l,2,...,Lx defines the length. 
The sum ^^j. j.,^ is taken only for nearest-neighbor (NN) 
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bonds. 

The values of t reported in the Hteratur^^HMIlII range 
from 2.4-2.7 eV for nanotubes, while i ~ 3 eV is typical 
in graphites. Although an accurate value of U is not yet 
known i n GN Rs, the estimate from polyacetylene, U ~ 
6-10 eV|22l22l might serve as a reasonable guess. Thus, 
we expect the ratio U /t to be of order one. 

Let us try to obtain the band structure of the armchair 
GNR within the tight-binding model first. For conve- 
nience, we consider the infinite length L^; ^ oo in the fol- 
lowing. Since the system is translational invariant along 
the x-direction, the hopping Hamiltonian simplifies after 
the partial Fourier transformation. 



c(r) 



,ik[x+Hy)] 



c(fc,y). 



(2) 



We omit the spin index a for a while. The shorthand 
notation y = (t;,A) is defined to label the sites within a 
unit cell. The geometric phases arisen from the underly- 
ing honeycomb structure are 



<5(y) = 2/3, 
Siy) = 0, 
'5(y) = 1/6, 
'5(y) = 1/2, 



for y=even and A = A, 
for y^even and K = B, 
for y=odd and h — A^ 
for y=odd and A = B. 



(3) 



The hopping Hamiltonian simplifies to decoupled two- leg 
ladders with finite length Ly labeled by momentum k in 
the a;-direction, 



(4) 



where ^{k) = [c{k, y; A), c{k, y; B)] is the fermion opera- 
tor for sublattices A and B and thus has 2Ly components. 
The reduced hopping Hamiltonian H{k) IS a '^Ly X '^Ly 
matrix. It is rather interesting that H(k) can be casted 
into supersymmetric (SUSY) form^^'^MSi 



H{k)^ 





Q 



(5) 



The submatrix that connects opposite sublattices is the 
supercharge operator. 



/t; ti 

ti tl ti 

ti t*2 h 

V 



(6) 



t*2 / 



where the complex hopping amplitudes are ti = — fe'*^/^ 
and t2 — — te*'^/^. One should not be surprised by the 
complex hopping amplitudes that come from the corre- 
sponding geometric phases kS{y). 

To solve the wave function, let us introduce a 2Ly- 
component spinor. 



'i>(y) = 



fsiy) 



(7) 



where fA/siy) are the wave functions on sublattices A/B. 
Making use of the SUSY form in Eq. ([5|, the coupled 
Harper equations are 



(8) 



It is known that the solution for E' ^ is supersymmetric 
and can be solved by combining the two Harper equations 
together Q^Q (fiA = E'^fA- Once the eigenstate lpa is 
obtained, the wave function on the other sublattice is 
= j^Q fA- Thus, the remaining task is to diagonalize 
the matrix Q^Q. Note that the above trick only works for 
E eigenstates while the E = states can be obtained 
by finding the null space of Q and alone. 

Before digging into details, we would like to make some 
general remarks. It is clear that, for each solution ipA, we 
can construct two wave functions on the other sublattice 
ifisby two choices of eigenenergies E = zt\E\. As a result , 
the energy spectrum is symmetric about E = and total 
wave functions of opposite energies E = ±\E\ only differ 
by an overall minus sign on one of the sublattices (but 
not on the other). 

With this general picture in mind, we now turn to the 
explicit form of the matrix Q^Q, that can be worked out 
by straightforward algebra 



/Vq~T2 Ti T2 

Ti Vo Ti 

T2 Ti Vo Ti T2 

T2 Ti Vo Ti T2 



V 







To Ti 



V0-T2J 

(9) 

where Vq = M^, Ti = 2t^ cos(fc/2) and T2 = t^. This ma- 
trix resembles the hopping Hamiltonian of a finite chain 
with the site potential Vg, nearest-neighbor hopping Ti 
and next-nearest-neighbor one T2. The eigenfunction sat- 
isfies 

T2[vA{y + 2) + VA{y - 2)] + Ti[fA{y + l) + <^^(y - l)] 
+Vo^a{v) = e\a{v), (10) 

where ?/ = 1, 2, iy with the appropriate boundary con- 
ditions 

(^^(0) = 0, ^^(L,y + l)-0, (11) 

V5a(-1) = -'/3a(1), VA{Ly + 2) = -ipA{Ly). (12) 

The first two boundary conditions arise from the open 
ends of the effective two-leg ladder and the last two con- 
straints comes from the change of the potential at the end 
sites. Note that the usual plane- wave solutions satisfy the 
bulk Harper equation in Eq. ( 10 ). Thus, we only need to 



form appropriate linear combination of these plane-wave 
solutions to match the boundary conditions. In the case 
we study here, the eigenstate is the simple combination 
of opposite momentum states with a relative minus sign, 
i.e. the wave function takes the usual sine form, 



(fiA{y) = sin((7,„y). 



(13) 
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where the magnitude of transverse momentum is quan- 
tized, Qm = rmr/(Ly + 1) with m = 1, 2, Ly. From the 
eigenstates, it is straightforward to compute the corre- 
sponding dispersions for each band, 

E = ±[Vo + 2TiCOsq,n + 2T2COs2q,n]^/^ 

= ±<[1 + 4cos(fc/2) cosg„i + 4cos^ (14) 

With the wave function on sublattice A and the energy 
dispersion, we can obtain the wave function on sublat- 
tice B by the supercharge operator as described before. 
However, a closer look would ensure us that the full wave 
function on the armchair GNR is far simpler than we ex- 
pected. 

The simplification arises from the fact that the super- 
charges Q and commute. As a result, the matrix 
QtQ = QQt share the same eigenstates as the matrices 
Q and . It is straightforward to show that the wave 
function (fAiv) is also an eigenstate of Q with complex 
eigenvalue. 



QifA = {t*2 +2tiCosqrn)ipA = \E\e' 



9(k; 



'Pa, 



(15) 



where the phase 6{k) = 6{k,q,n) of the complex eigen- 
value is 



„j0(k) _ *2 + 2^1 cos g„ 



(16) 



Therefore, the wave function on the sublattice _B is a 
duplicate of with a momentum-dependent phase shift. 



<f,iy) = ±e''^'^^ sin(g„y). 



(17) 



The ± signs arise from the signs of the energy, corre- 
sponding to antibonding and bonding bands respectively. 
Finally, after attaching appropriate renormalization fac- 
tor, the full eigenstate wave function is labeled by the 
quantum number k = (fc, g^js), including the longitu- 
dinal momentum k, the transverse momentum qm and 
antibonding/bonding index s = ±1. The explicit form of 
the wave function is 



'i>k(y) 



1 



\/Ey 



sin{qmy) 
ggi6»(k) sm{q^y), 



(18) 



Note that the dependence of the longitudinal momentum 
k only appears through the relative phase 0(k) between 
wave functions on two sublattices. This simple analytical 
form of the eigenstates allows us to map the armchair 
GNR to effective theory in the low-energy limit and study 
the correlation effects analytically. 

For the width of odd Ly, the transverse momentum 
goes through the particular value q^ — 7r/2, rendering 
the energy dispersion completely flat at energy E — ±t 
in the whole Brillouin zone. To obtain the wave function, 
we only need to evaluate the phase shift, e'^^± — m + 
2ti cos(7r/2)]/t = -e"''^/^. Therefore, the wave functions 
for the flat bands E — ±t are 



*.±(y) 



1 



sin(7r?//2) 
^g-»fc/3sin(7r?//2) 



(19) 



Since all states with different momentum k are exactly 
degenerate, it is possible to construct the local Wannier 
orbital with the same energy, 



1 



1 



sin(7r2//2) 
=Fsin(7ry/2) 



(20) 



Note that the momentum-dependent phase shift 9p± = 
—k/3 cancels the relative geometric phase k[S{y,B) — 
6{y, A)] leading to extremely simple local Wannier orbital 
at X = xo- Repeatedly applying the lattice displacement 
operator on one Wannier orbital, all orbitals at dif- 
ferent locations can be constructed. Since [Tj;,Ht] ~ 0, 
all the orbitals share the same energy and form a flat 
band. This is the one-dimensional analog of the Landau 
level degeneracy for two-dimensional electrons in mag- 
netic field. The peculiar edge topology at nanoscale re- 
places the role of the magnetic field in 2D and quenches 
the kinetic energy of the carriers. 

At first sight, it is rather counterintuitive that the local 
Wannier orbital cannot move around by quantum hop- 
ping. The static nature is due to perfect destructive 
quantum interferences which make hopping amplitudes 
from different sites cancel each other. Thus, the open 
boundaries of armchair shape play a crucial role for the 
birth of the Wannier orbitals. Furthermore, since the 
wave function is not zero only at odd y coordinates (see 
Fig. [T|), it is clear that the adjacent orbitals have zero 
overlap. Thus, the Mielke-Tasaki mechanism does not 
work to couple neighboring orbitals magnetically. 

Let us concentrate on the flat-band regimes E = ±t for 
the armchair GNR with odd Ly. Due to the particle-hole 
symmetry for the Hubbard model, the low-energy physics 
is dictated by the doping level disregarding whether it is 
electron or hole doped. Thus, it is convenient to intro- 
duce the positive-definite doping level Xd = \{n) — 1|, 
where (n) is the average electron number per site. The 
lower and upper bounds of the doping rate Xd for the 
fiat-band regime are obtained from the Fermi momenta 
km of the dispersive bands intersecting the fiat band. 
For those dispersive bands, the Fermi momentum sat- 
isfies cos(/cm/2) -I- cos(gm) = 0, which leads to fc,„ — 
27r — 2q„i, where m — Ly,Ly — l,...,{Ly + 3)/2: there 
are {Ly — l)/2 pairs of Fermi points crossing the fiat 
band. Thus, the system is in the fiat-band regime for 
a;<i,min < Xd < a;d,max, where 



1 



XdA 



Xdj 



— T' 

T:Ly ^ 



XdA 



1 



4L,' 



1 



(21) 



1 

4 

_ 1 3 

Ly " 4 iLy 

Therefore, the fiat-band regime shrinks as the transverse 
width increases and eventually goes to zero in the two- 
dimansional limit. This trend highlights the importance 
of finite transverse width of the system and why the fiat- 
band physics is no longer the dominant player in 2D. In 
the following, we will try to write down the effective field 
theory for both the fiat and dispersive bands in weak 
coupling. 
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III. WEAK-COUPLING ANALYSIS 

Even though wc have derived the analytical form of 
wave functions in armchair GNR, it is still quite compli- 
cated to write down the effective field theory. In the flat- 
band regime, after integrating out gapped bands, there 
remains a flat band with intersecting dispersive bands. 
Note that the lattice fermion can be decomposed into 
eigenstates. 



Ca{x,y,h) 



1 

/ L„ 



ikx \ ^ 



(y) (fc), (22) 



where the extra geometric factor is included in the mod- 
ified wave function (jimsiy) = e^'^^^^^^msiy)- In the 
low-energy limit, one can approximate all intersecting 
bands with linear dispersions, the above expansion is 
then greatly simplified. 

Let us work out the chiral-field expansion in the Hat- 
band regime at i? = i as an example. In that case, the 
effective low-energy theory is described by the flat band 
and the intersecting dispersive bands of the antibonding 
sector s = 1. Thus, the lattice fermion can be decom- 
posed into the flat-band and pairs of chiral field opera- 
tors. 



Ca(r) ~ <?!)p(y)V'F„(a;) -f 



(23) 



where P — R/L — +/— represents the chirality. The 
modified (including the geometric phases) wave function 
for the flat-band orbital is 



1 



sin(7r?//2) 
— sin(7rj//2) 



(24) 



and those for the right/left-moving plane waves at Fermi 
point ±fcm are 



0p™(y) 



1 



gip»c„[2/3 + .Blsin(g„y) 



(25) 



A/l 



where the shorthand notation is introduced S 
5{y, A/ B). Substituting the chiral-field expansion into 
the lattice Hamiltonian, one can easily derive the effec- 
tive field theory in the flat-band regimes. 

Since the density of states is divergent for the flat band, 
the dispersive bands can be dropped to the leading order 
approximation. The effective Hamiltonian, keeping the 
flat band only, is rather simple. 



Hp ~ L/p ^ np^(a;)np^(a;), 



(26) 



where ripA^) is the number density of each spin flavor 
and Up is the effective on-site interaction for the flat-band 
orbitals, which can be computed by projection onto the 
flat band, 



Up = U 



y 



u 



(27) 



Note that the kinetic energy is quenched in the flat band 
and thus the ground state contains no quantum fluctua- 
tions. To avoid the cost of Up, the ground state consists 
of the maximum number of singly-occupied Wannier or- 
bitals, which leads to local magnetic moments. Since 
there is no overlap between adjacent orbitals, these mag- 
netic moments are free and give rise to a large ground- 
state degeneracy. 

To lift the large degeneracy of the ground state, inter- 
action with the itinerant carriers in the dispersive bands 
must be included. While the effective Hamiltonian for 
the interaction contains other terms, the terms to play a 
key role are the exchange interactions which couple the 
local moments and the itinerant carriers. 



(28) 



where Sj7(x) is the spin density operator for the local 
moments and Smix) = SH,„(a;) + Si„(x) is the spin den- 
sity of itinerant carriers in the crossing band m. The 
exchange integral is given by. 



Jp„ = 2^0;(y)<^,(y')0P 



.(y)0;™(y')^y,y'- (29) 



Using Eqs. ( 24 1 and ( 25 ) , we obtain the exchange integral 
due to the on-site interaction Vy^yi — Sy^y'U; it has a 
rather simple form. 



y 

2U 



2sm\q^y) = (30) 



(i„ + 1)2 

^ y ' y=odd 

The subscript P — R/L is dropped because the couphng 
strength does not depend on the chirality. Thus, the on- 
site interaction induces a ferromagnetic coupling between 
the local moments in the flat band and the itinerant spins 
in the dispersive bands. 

The exchange coupling Jm tends to align the local mo- 
ments from the flat band because it does not cause any 
extra kinetic energy. The ferromagnetically aligned mo- 
ments act back on the itinerant carriers and induce fi- 
nite polarization in the dispersive bands. The interacting 
Hamiltonian Hp + H , therefore shows interesting two- 
step flat-band ferromagnetism - electronic correlations 
in the flat band give rise to local moments without direct 
exchange coupling, while the presence of gapless itiner- 
ant carriers mediates the ferromagnetic order. The sig- 
niflcant feature of the armchair GNR is that, even within 
the one-orbital Hubbard model without any magnetic im- 
purity nor additional localized levels, the electronic cor- 
relations give rise to both local moments and itinerant 
carriers simultaneously due to the peculiar topology of 
the open edges. 

It is also interesting to consider the flnite-size effect 
arisen from the length of the armchair GNR. If one 
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FIG. 2: (Color online) Doping dependence of the energy dif- 
ference A{xd,SS) between higher-spin and lowest-spin states 
in the armchair GNR with Lx = 2 and Ly = 5 for U/t = 4. 
The squares, triangles, and diamonds represent the data for 
dS = S — So = 1, 2, and 3, respectively. 



TABLE I: Energy difference per unit cell A(a;d) between the 
ferromagnetic ground state and the paramagnetic one in the 



armchair GNR with Lx 



2 and L„ 



5 for the flat-band 



regime. Here S is the total spin of the ground state and the 
hopping amplitude is chosen t = 3 eV. 





5" 


U/t = 2 


U/t = 4 


U/t = 8 


0.25 


3/2 


-44 meV 


-54 meV 


-57 meV 


0.30 


2 


-115 meV 


-111 meV 


-87 meV 


0.35 


3/2 


-48 meV 


-52 meV 


-34 meV 



imposes the periodic boundary conditions along the x- 
direction, the system becomes a short segment of arm- 
chair nanotube. In this case, only when the quantized 
momenta — 21ti/Lx {I = 0,1,..., — 1) coincide 
with the Fermi points ±fcm, the gapless itinerant car- 
riers are present and the ferromagnetic ground state 
is realized. For open boundary conditions with finite 
Lx, the situation is much more complicated; is no 
longer good quantum number and the band structure 
can be deformed by finite L^. However, from numeri- 
cal calculations for the tight-binding Hamiltonian Ht , we 
have found that the energy spectrum around the flat- 
band level £' = ±t is not affected by finite and 
can be accurately approximated by the quantization rule 
kx = Itt/{Lx + 1) {I — 1,2, ...,Lx). When the quantized 
momenta k^ coincide with the Fermi points km, ferro- 
magnetism sets in with the help of these gapless itinerant 
carriers. Therefore, depending on specific choice of L^, 
the ground state of the armchair GNR in the flat-band 
regime can be ferromagnetic or Curie-like paramagnetic. 




10 
(c) site index /' 



FIG. 3: (Color online) Spin polarization profile {S^(r)) of 
the ferromagnetic ground state for the armchair GNR with 
Lx ~ 2 and Ly = 5. It is at the optimal doping x'^ — 0.30 
with the interaction strength (a) U/t = 4. (b) U/t = 0+. The 
values of {S^{r)) are represented by the areas of the shaded 
circles. It is remarkable that the profiles with intermediate 
and weak interaction strength are almost identical, (c) The 
values of {S'^(r)) at each site for f//t = 0^, 2, 4, 8 represented 
by solid circles, squares, triangles and diamonds. 



IV. NON-ABELIAN DENSITY MATRIX 
RENORMALIZATION GROUP 

To check the validity of the weak-coupling scenario and 
see whether it survives for the realistic coupling regime, 
we choose the non- Abelian DMRG method.^'^ It IS im- 
portant to emphasize that, to look for the higher-spin 
ground state, the non- Abelian approach is more powerful 
and convenient compared to the conventional DMRG'^^ ^ 
because the former makes the full use of the spin SU(2) 
symmetry. Employing the non- Abelian DMRG method, 
we can compute the energy difference between the higher- 
spin (ferromagnetic) state and the lowest-spin (paramag- 
netic) state. 



A{xd,SS) = Eo{xd,S) ~ EQ{xd,So), 



(31) 



where Eo(xd, S) is the lowest energy in the subspace with 
the doping rate Xd and the total spin S. Furthermore, 
SS = S — So, where 5*0 — 0, 1/2 denotes the lowest spin 
depending on whether the number of carriers is even or 
odd. The calculation is performed for the system with 
Ly — 5 and L^ — 2, for which the quantized momenta 
kx = It:/ 3 coincide with the Fermi points of the dispersive 
bands and therefore the weak-coupling theory predicts 
the ferromagnetic ground state. The number of SU(2) 
multiplets kept is up to 450, typically corresponding to 
1000-3000 U(l) states. The truncation error is of order 
10^^ or less, and the results are extrapolated to the limit 
of zero truncation error. 

Figure [2] shows the energy difference A{xd,SS) as a 
function of the doping rate Xd and the spin SS for 
U/t — 4. We do find numerically the ferromagnetic 
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(higher-spin) ground state in the flat-band regime, Xd = 
0.25, 0.30, 0.SS.-^The results with U/t = 2, 8 (not shown 
here) also show a similar doping-rate dependence, sup- 
porting the ferromagnetism for the flat-band regime. 

Collecting all data for A{xd, SS) together, we can de- 
termine the ground state at each doping level and its en- 
ergy gain per unit cell, A^x^), which is the energy differ- 
ence between the higher-spin ferromagnetic ground state 
and the paramagnetic one with lowest spin, 

A{xd) = -^[Eoixd) - Eoixd^So)], (32) 

where Eo^Xd) is the ground state energy at doping level 
Xd- The results for the flat-band regime (0.2 < Xd < 0.4 
for Ly — 5) and U/t — 2, A, 8 are summarized in Table |l] 
The optimal doping occurs at x*^ = 0.3 as predicted by 
the weak-coupling theory. Thus, the non-Abelian DMRG 
results support the carrier-mediated ferromagnetism pre- 
dicted from the analytical approach in weak coupling 
even when the interaction strength is in the intermediate 
regime. 

To further verify the role of the flat band, we also calcu- 
late the profile of spin density in the ground state. Since 
the interaction strength is now larger than the hopping 
amplitude, one may guess any peculiar feature in the 
band structure should be suppressed. Figure |3] shows 
the results at the optimal doping Xd — 0.3 with total 
spin S = 2. Remarkably, the spin polarization for finite 
U/t = A has a similar profile to that in the weak-coupling 
limit U 0^ obtained from the cigcn-wavefunctions of 
the tight-binding model Ht- The result clearly indicates 
that the flat-band orbit als still play a signiflcant role 
in the ferromagnetic ground state for finite U/t. The 
physical picture developed in weak coupling thus applies 
rather well and extends smoothly to the intermediate- 
and strong-coupling regime. 

As the momentum kx is discretized in the system with 
finite Lx, it is important to see how the properties of 
the system change depending on L^- For armchair GNR, 
the large unit cell and peculiar nature of the fiat-band 
states lead to slow convergence of the DMRG calculation 
and make it difficult to treat the system with larger L^, 
unfortunately.'^ Nevertheless, we have performed numer- 
ical calculation for other one-dimensional-lattice models 
which have essentially the same band structure as that 
of armchair GNR. We have then found that the models 
with larger number of unit cells indeed exhibits the itin- 
erant ferromagnetism with the properties expected from 
the weak-coupling theory in Sec. even including the 
peculiar finite-size effect. The results will be reported 
elsewhere.l^ Furthermore, we also emphasize that, as 
we will see in Sec. |V] the first-principles calculation for 
infinite-length armchair GNR also shows the itinerant 
ferromagnetism for the flat-band regime, in accordance 
with the weak-coupling prediction. These observations 
support that the carrier-mediated ferromagnetism found 
in this section would survive for larger and connect 
to the thermodynamic limit oo. 




FIG. 4: (a) Band structure of the undoped, infinitely long 
armchair GNR with Ly = 5. (b) Doping dependence of the 
spin polarization per unit cell, (c) Energy difference per unit 
cell A{xh) between ferromagnetic and paramagnetic states as 
a function of the hole doping. 



V. LOCAL SPIN DENSITY APPROXIMATION 

The complimentary approaches of the weak-coupling 
analysis and the non-Abelian DMRG method establish 
the ferromagnetic ground state in armchair GNR within 
the Hubbard-type model. To treat GNRs, however, one 
must take account of the realistic band structure be- 
yond the tight-binding approximation as well as the effect 
of unscreened Coulomb interactions. For the purpose, 
we have carried out first-principles calculation of the lo- 
cal spin-density approximation within density functional 
theory. The self-consistent band structure calculations 
under lattice optimization were performed usingthc full- 
potential projected augmented wave methoc f " ' as im- 
plemented in the VASP package.!^^!^ 

Figure Qa) shows our result of the band structure of 
the armchair GNR with Ly — 5 at half filling. The nu- 
merics show that the band structure of the undoped arm- 
chair GNR is more or less similar to the tight-binding re- 
sults. There exists a narrow band with bandwidth ~ 0.4 
eV located at 2.5 ~ 2.9 eV below the Fermi level crossed 
by two itinerant bands. Furthermore, one can find sim- 
ilar structure in the unoccupied counterpart though the 
particle-hole symmetry does not hold for this case. 

Although the lower "fiat band" has a finite band- 
width 0.4 eV) and thus is not fiat anymore, it is 
still massive enough compared with the other two dis- 
persive bands. As a result, the two-step carrier-mediated 
ferromagnetism still works as will be described below. 
Figure presents the magnetic moment per unit cell 
at different hole doping levels Xh = I — {n). Upon hole 
doping, the magnetic moment goes up and reaches its 
maximum, ~ l.l/is per unit cell, at the optimal hole 
doping a;^ ~ 0.28. In the mean time, as shown in Fig. 
|4[c), the energy gain of the ferromagnetic state per unit 
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romagnetic regime at finite interaction strength U. Con- 
sidering the GNR sHghtly outside the flat-band regime, 
although costing higher kinetic energy, it is still prefer- 
ential to fill in particles/holes in the flat band since the 
exchange energy is lowered. Thus, we expect that the fer- 
romagnetic ground state can exist even outside the flat- 
band regime. The variational calculations indeed show 
that the above understanding is correct and the ferro- 
magnetic phase exists in a wider range than the flat-band 
regime. Therefore, it is expected that the inclusion of 
the realistic Coulomb interaction will give rise to similar 
enhancement of the flat-band ferromagnetism as demon- 
strated in Fig. |4]here. 



FIG. 5: (Color online) (a) Spin decomposed band structure 
of the infinitely long GNR with = 5 at the optimal hole 
doping xj^ ~ 0.28. (b) Top and (c) side views of the spin 
density distribution (shown in grey color). The yellow (light) 
and red (dark) spheres denote C and H atoms respectively. 



cell raises noticeably, especially near the optimal doping, 
up to ~ 37 meV/cell. These results clearly support the 
emergence of the ferromagnetism. 

The obtained value of the magnetic moment ~ 1.1 /is 
per unit cell indicates that the nearly-flat band is al- 
most fully-polarized at the optimal doping. This feature 
is also shown by the spin decomposed band structure in 
Fig. ^a.). The bandwidth of the nearly-flat band at 
the optimal doping is slightly suppressed to ^ 0.3 eV 
because of the reduced density. The exchange splitting 
between opposite spins is about 0.5 eV. Figures [sjb) and 
(c) demonstrate respectively the top and side views of 
the spin density distribution at the optimal doping. It 
is truly remarkable that the proflle of the spin density 
is almost identical to the Wannier orbital in the weak- 
coupling limit (see Fig. [T]). The nodes caused by de- 
structive quantum interferences can be seen clearly in 
the numerics. We emphasize that not only the dumb- 
bell shape of the spin density emerges as predicted, the 
optimal doping concentration and the size of the spin po- 
larization realized also agree with the prediction of the 
weak-coupling theory. 

As the hole doping increases further, the magnetic mo- 
ment as well as the energy gain of the ferromagnetic state 
are suppressed. The ground state eventually turns non- 
magnetic at larger hole doping. It is interesting to notice 
that the ferromagnetic ground state, as shown in Fig. |4] 
exists for 0.15 < < 0.42, which is wider than the 
flat-band regime between Xmin = 0.2 and Xmax = 0.4 
predicted by weak-coupling theory for Ly — 5 armchair 
GNR. The ferromagnetic phase obtained by the local 
spin-density approximation then seems to persist even 
slightly outside the flat-band regime. This shall not be 
surprising since the interaction strength is no longer weak 
in comparison with the kinetic energy. To address this 
issue, one can calculate the variational energy of the Hub- 
bard model in the armchair GNR^ to pin down the fer- 



VI. DISCUSSIONS AND CONCLUSIONS 

Here we would like to elaborate on several subtle points 
which were not discussed in previous sections. First, we 
note that the weak-coupling theory in Sec. |III| developed 
for the Hubbard model in the armchair GNR is pretty 
robust against lattice distortions. In practical graphite 
network, the hopping along horizontal and tilted vertical 
bonds is expected to be slightly different. Though the 
deviation shifts the flat band from E = ±.t, the bands re- 
main flat and the same mechanism of the ferromagnetism 
applies. 

Furthermore, the profile of the short-range interaction 
does not seem to do much harm either. Following steps 
similar to those in Sec. |III| we can compute the exchange 
coupling due to nearest-neighbor interaction V± (tilted 
vertical bonds) and V|| (horizontal bonds). Since the 
product of the flat-band wave function (j)p{y)cj)p{y') = 
if the coordinates are connected by a tilted vertical bond, 
V± does not give rise to any exchange coupling. On the 
other hand, the nearest-neighbor interaction along the 
horizontal bonds V|[ gives rise to non-vanishing exchange 
coupling. 



Jrr 



y 

2Vn 



{Ly + 1) 

Ml 



Ly+l 



— (-2cosfc„j) ^ sw?{qmy) 

y—odd 

COS km- (33) 



Here we have used the expression for the interaction 
Vy^y' ~ <5y,y'^| with the notation y = {y,B/A) for 
y — {y,A/B). As expected, the exchange couplings 
arisen from right-/left-moving fields in the same disper- 
sive band are the same and thus the subscript P = R/L 
is dropped. We thereby find that, in the presence of 
the nearest-neighbor interaction Vj_ and V|[ , the exchange 
coupling in Eq. ( 28 1 becomes 



1 



L,, 



(t/ - V|| cos km)- 



(34) 
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Since , V|| < U is often expected, the exchange cou- 
pling is still ferromagnetic and the picture does not 
change. 

The above calculations can be generalized to 
the screened short-ranged interaction. Suppose the 
spatial profil e of t he screened interaction goes as 
exp(— a;/^)/-\/a;^ -I- Iq, where Iq is a short-range cutoff 
(comparable to the lattice constant) and ^ is the length 
scale of the short-ranged interaction. Ignoring the com- 
plicated form factor due to detail orbital overlapping, 
the exchange coupling takes the general form Jm{x) ^ 
exp{—x/(,)cos{kjnx)/x. As expected, the exchange cou- 
pling decreases as the distance is far apart. Furthermore, 
the oscillatory factor cos{kmx) makes the couplings to the 
dispersive bands with different signs and tends to cancel 
each other. It will further suppress the effects of the 
exchange coupling beyond the nearest neighbors. This 
trend is in agreement with our first-principles calcula- 
tions where the true long-ranged Coulomb interaction is 
included. 

To realize the fiat-band ferromagnetism in armchair 
GNR, the crucial challenge lies in how to achieve the ap- 
propriate doping level. One of graphene's superior prop- 
erties is its pronounced ambipolar electric field effect .EEEl 
By applying gate voltages, the charge carriers can be 
tuned between electrons and holes with concentration up 
to 10^^ cm~^. Even so, it is unlikely that the exter- 



nal gate voltage alone can pour enough electrons/holes 
into the system to reach the flat-band regime. An- 
other route to dope GNR is via chemical doping. It was 
demonstratecPthat the chemical dopants in the substrate 
can markably change the carrier density. Perhaps the 
combination of both methods can be even more efficient. 

In conclusion, by combining the weak-coupling anal- 
ysis, the non-Abelian DMRG technique, and the first- 
principles calculations, we show how ferromagnetism oc- 
curs in armchair GNR - electronic correlations give rise 
to magnetic moments in the fiat band and the itinerant 
carriers in the dispersive bands mediate ferromagnetic 
coupling between these uncoupled moments. Recently, 
there are proposald^^'^ to use GNR to build transistors 
and spin qubits. While these proposals take care of many 
realistic issues, the electronic correlations are ignored. 
Our study here show that electronic correlations in GNR 
can bring up surprises such as the carrier-mediated fiat- 
band ferromagnetism. Therefore, it is crucially impor- 
tant to include the correlation effects when we try to 
realize these proposals into devices. 
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